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Abstract 

We assume that, in equilibrium, nuclear matter at reduced density and moderate finite temper- 
ature, breaks up into many fragments. A strong support to this assumption is provided by data 
accumulated from intermediate energy heavy ion collisions. The break-up of hot and expanded nu- 
clear matter according to rules of equilibrium statistical mechanics is the multifragmentation model. 
The model gives a first-order phase transition. This is studied in detail here. Phase-equilibrium 
lines for different degrees of asymmetry are computed. 

PACS numbers: 25.70Mn, 25.70Pq 
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I. INTRODUCTION 



Nuclear matter is a hypothetical very large system of nucleons where the Coulomb effects 
of protons are switched off. Such a system is expected to have features of hquid-gas phase 
transition. We consider here the equation of state of symmetric and asymmetric nuclear 
matter at temperature between 4 and 10 MeV and at less than half of normal nuclear 
density. We assume that at equilibrium at finite temperature (three to tens of MeV) and 
low average density, nuclear matter breaks up into fragments, each with normal nuclear 
density. Strong support to this assumption comes from data on heavy ion collisions but it is 
also supported by theoretical modelling. For example it can be easily shown (see section IV 
in using Skyrme type interaction, that the free energy of uniformly stretched nuclear 
matter is very significantly lowered if the matter is allowed to split into many fragments, 
each with normal nuclear density. This is the multifragmentation model. We use this model 
to study thermodynamic properties of nuclear matter, particularly phase-equilibrium lines 
(the lines of co-existence of liquid and gas phases) in the p — T plane for both symmetric 
and asymmetric matter. 

This is an extension of the model described in our earlier work where only one kind of 
particles was considered. This one kind of particles however formed clusters whose properties 
were patterned after actual finite nuclei. While we hope that the present article is self- 
contained we will refer to this earlier work for elucidation of some points. There is a large 
number of publications on equation of state and phase transitions in nuclear matter. Ref. 
[sl comes closest to the spirit of this paper. While there are quite a few common features 
with this work there are also some differences and we highlight some other aspects. Phase 
transition in nuclear matter using mean-field theory was studied over many years and we 
can not attempt an adequate bibliography here. We mention two papers which critically 
looked at asymmetric nuclear matter and received a great deal of attention in very recent 
times 0, [sl. Both of these use mean- field theories and overcome the difficulty of instability 
through Maxwell construction. The multifragmentation approach is very different. It is 
more directly related to actual observables but in its present form it can only be trusted in 
a low density regime. But there is no need for Maxwell construction. 
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II. THE FORMULAE 



We briefly review the grand canonical model for multifragmentation 6|. Let the numbers 
of neutrons and protons in the dissociating system be A^o and Zq respectively. At finite 
temperature and in subnormal densities, these will break up into all possible composites 
each with some neutrons(A^) and protons (Z)(mass number A = N + Z). We always use 
the subscripted A^^o, Zq to refer to the very large system whose thermodynamic properties 
are being investigated whereas A^, Z refer to composites which can be small or large. The 
properties of the composites are determined by the basic two-body interactions These prop- 
erties are utilized in the model but interactions between composites are neglected (except 
through excluded volume effect; see discussion later) by appeahng to the short range nature 
of nuclear forces. This limits the validity of the model to low densities. Here we will restrict 
ou.. ...estigafou to de,.,t.s ,/p„ to 0.5 o. less whe.e p. . the ,>o™a. ..c.ea.. density. This 
IS the customary practice [3J. 

If the neutron chemical potential is fin and the proton chemical potential is fip, then 
statistical equilibrium implies that the chemical potential of a composite with A^ neutrons 
and Z protons is /i„A^ + fipZ. The following are the relevant equations for us. The average 
number of composites with A^ neutrons and Z is (/3 = 1/T) 

M = e^>^-'^^^^^^u;j,,z (1) 

Here uj^^z is a one body partition function for the composite (A^, Z). It is a product of two 
factors; one arising from the translational motion of the composite and another from the 
intrinsic partition function of the composite: 

u;n,z = ^(27rmT)3/2A3/2 x ZN,z{int) (2) 

Here Vj is the volume available for translational motion; Vf will be less than V, the volume 
to which the system has expanded at break up (excluded volume correction). We use 
Vf = V — Vq , where Vq is the normal volume of nucleus with Zq protons and A^o neutrons. 
The quantity ZN^zi'int) depends upon the intrinsic properties of the composites and contains 
all the nuclear physics. 

We list now the properties of the composites used in this work. The proton and 
the neutron are fundamental building blocks thus Zi^Q{int) = ZQ^i{int) = 2 where 2 



3 



takes care of the spin degeneracy. For deuteron, triton, ^He and ^He we use Zij{int) = 
{2sij + 1) exp{—Peij{gr)) where eij{gr) is the ground state energy of the composite and 
(2sjj + 1) is the experimental spin degeneracy of the ground state. Because we are modehng 
a system where protons do not carry any charges the ground state energy of ^He is taken to 
be that of the triton and the Coulomb energy is subtracted from the experimental energy 
of the alpha particle. These modifications make insignificant changes. Excited states for 
these very low mass nuclei are not included. For mass number a = 5 and greater we use the 
liquid-drop formula. This reads 

Zij{int) = exp[ — ^] (3) 

Here Fij is the internal free energy of species 

F„. = -Woa + a(T)a2/3 + s^^^^^ - (4) 

a eo 

The expression includes the volume energy, the temperature dependent surface energy and 
the symmetry energy. The values of the paprameters are taken from fsl . The term rep- 
resents contribution from excited states since the composites are at a non-zero temperature. 
For nuclei with A=5 we include Z=2 and 3 and for A=6 we include Z=2,3 and 4. For 
higher masses we compute the drip lines using the liquid-drop formula above and include all 
isotopes within these boundaries. 

There are two equations which determine /i„ and fip. 

No = J2 Ne^^-''^^^^^UN,z (5) 

Zo = E^e^""''+^"^^^iv,z (6) 

We want to point out the following feature of the grand canonical model. In all ujn,z^^ in 
the sum in the above two equations, there is one common value for Vf (see eq.(2)). We 
really solve for Nq/Vj and Zo/Vf. The values of fin or ftp will not change if we, say, double 
Nq,Zo and Vf simultaneously provided the number of terms in the sum is unaltered. We 
then might as well say that when we are solving the grand canonical equation we are really 
solving for an infinite system (because we know that fluctuations will become unimportant) 
but this infinite system can break up into only certain kinds of species as are included in the 
above two equations. Which composites are included in the sum is an important physical 
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ingredient in the model but intensive quantities like /3, /i depend not on Nq, Zq but on No/Vf 
and Zo/Vf. 

The choice of which nuclei are included in the sum of the right hand side of eqs. (5) and 
(6) needs further elucidation. We can look upon the sum on and Z as a sum over A and a 
sum over Z. In principle A goes from 1 to oo and for a given A, Z can go from to A. Here for 
a given A we restrict Z by the drip lines. Comparisons with calculations where restrictions 
by drip lines are not imposed (as in the Copenhagen statistical multifragmentation model) 
showed that restrictions by driplines generate imperceptible differences jo]. De et al [st reach 
similar conclusion. Let us now consider the restriction on A. In principle this should be oo 
but for practical calculations one needs to restrict this to a maximum value that we label as 
Amax- Earlier calculations with one kind of particles showed that with A^ax = 200 features 
of liquid-gas phase transition are not revealed (see Fig. 14 in but a high value of A^ax 
at 2000 produces a nearly a perfect model of phase transition (elaborated in much larger 
detail in Q and in Q). 

III. SIGNATURES OF PHASE TRANSITION IN THE MODEL 

We now demonstrate that the multifragmentation model predicts first order phase tran- 
sition. There are three signatures we will dwell on. Pressure in the model is given by 
p = T ^y^'^ = T = Tpj ^"^ . We plot results as function of p rather than pf the 
connection being pf = We have p = Pn + Pp- We need an asymmetry parameter. We 

use both No/Zq and u = 

We show in Fig.l p — p curves for Nq/Zq = 1.4 where the values of A^ax are 200, 400, 
600, 800 and 1000. The temperature used is T = 6.5 MeV. For all five choices of A^ax 
pressure against p initially rises quite sharply and then flattens out considerably. The 
initial stage of fast rise of pressure with density is the gas phase. Here the results do not 
matter whether A^ax is 200, 400 or larger. The reason will become clearer later (it is 
explained in detail in j^). The flattening which follows depend on A^ax but above a large 
enourfi value of A^ax will not change. For one kind of particles this is reached around 
2000 [21]. However, the choice of Amax=QOO is good enough for at least a semi- quantitative 
estimate of various thermodynamic properties of nuclear matter and we will present results 
for this value although we did some calculations with other choices of Amax also. The 
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flattening happens slightly beyond p/po=0.1. We show results up to p/po = 0.5 arguing 
that the excluded volume correction for interactions between composites becomes worse 
with increasing density. 

The rise of pressure at small density followed by a flattening of p with increasing density 
is a signature of first-order liquid-gas transition. We have shown results for T=6.5 MeV. 
Beyond a certain temperature the flatness will disappear showing that there is no more 
phase transition in the domain p/po < 0.5. Similarly the flattening of p disappears beyond 
some value of Nq/Zq. The liquid-drop parameters we are using give us for large nuclei the 
drip-line at N/Z (and of course Z/N) about 2. Hence for larger values of Nq/Zq the system 
can not stay together even at T=0. Then we will have a system which has a bound core 
but always many free nucleons which will dominate the thermodynamic properties of the 
system. This is not a system we want to study. Hence in this work we constrained ourselves 
to system whose Nq/Zq spans 1.0 to 1.8. The upper limit is indeed a highly asymmetric 
system. 

Below the density where phase transition sets in, the system is in pure gas phase. At 
phase transition point some liquid will be formed and the fraction of nucleons in the liquid 
phase will grow at the expense of the gas particles as the density increases. This can actually 
be followed. One also gets a functional definition of what constitutes the gas particles. Here 



our identification is very different from what is concluded in [3[ but very similar to what is 
found in our earlier work with one kind of particles 2|. 

Lastly, in one component model there is just one p which stays constant throughout the 
co-existence region. Now there are two chemical potentials pn and pp. How do they behave? 



IV. WHAT CONSTITUTES THE GAS AND WHAT CONSTITUTES THE LIQ- 
UID? 

The quantity < ua >= Y.n+z=a < fT'N,z > is the average number of composites with 
mass number A. The quantity A < ua > / Aq gives the fraction of particles tied up in com- 
posites with mass number A. This is plotted in Figs. 2 and 3 for Nq/Zq=\.^ and Nq/Zq=\.'^ 
respectively. First concentrating on Fig.2 (T=6.5 MeV) we see that at density p/po=0.1 the 
nucleons are bound in composites <50. These particles constitute the gas phase. At density 
p/po =0.3 some heavy composites with A ^ Amax begin to form and the probability of such 
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heavy particles (with A between Amax and Amax-^00) begins to increase (at the expense of 
the hght particles) as the density increases. This is a clear evidence of co-existence. We thus 
consider light particles {A < 70) to be gas and heavier particles (with A between A^ax and 
^max-135) to be liquids. Fig.3 displays similar physics but for Nq/Zq=1.8 : all gas particles 
at p/po=0.1 and mixture of gas and liquid at p/po=0.38. 

We note that even the gas phase in fragmentation model is quite complicated. It is not 
just neutrons and protons but other light nuclei as well. In addition, during co-existence the 
isotopic content of the gas phase changes continuously as the volume of the container, i.e. 
density p/po changes. This is called isospin fractionation and is well-known in literature. 
We will briefly come back to this aspect later. 

V. CHEMICAL POTENTIALS 

In numerical work involving one kind of particles only it was demonstrated that in 
the limit of Amax — * oo a constant value of p will be achieved in the co-existence region. 
This value could be obtained by extrapolation. In the present case there are two chemical 
potentials. For Nq/Zq 7^ 1, p„ 7^ pp. For N/Z=1A and temperature 6.5 MeV we show in 
Fig. 4 the evolution of p„ and Pp as a function of density. One notices that both p„ and Pp 
change rapidly in the gas phase and then tend to a constant value. In the limit A^ax — > cxd we 
expect they will become constants. We also plot in the same figure p = jj^^fz^f^n + Ivq+zJ^^p- 
The p so defined has a meaning at the three limits:-!, and +1 for asymmetry parameter 
^ = nI+zI interesting to note that p tends to constant value faster than either p„ 

or Pp. 

VI. CO-EXISTENCE LINES 

Figure 5 shows that as the temperature increases phase co-existence finally disappears 
(from the region p/po ^ 0.5). We have shown this for No/Zq=1.0 but this is also true for 
asymmetric systems provided the asymmetry is not too large as explained earlier. 

We show in Fig. 6 the p — p curves at T=6.5 MeV for three systems with Nq/Zq=1, 1.4 
and 1.8. We can identify from the figure points A, B and C on these curves where co- 
existence sets in. The values of pressure at these points give us p values for co-existence 
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at this temperature for these Nq/Zq values.. This is not strictly correct. The values of 
p increase slightly as one moves towards higher density. This is because with Amax=Q^^ 
we have not reached asymptotic limits yet. However, this is adequate for our purpose. 
Repeating this analysis for different temperatures we get co-existence lines m T — p plane 
for nuclear matter with different asymmetries (Fig. 7). Notice that while the co-existence 
lines for differing asymmetries are different they are quite similar. As usual, points to the 
left and above the co-existence lines are in the liquid phase and points to the right and below 
are in the gas phase. 

The highest point of a co-existence line in the T — p plane usually identifies critical values 
Tc-iPc [3- This is not true in Fig. 7. As we consider higher temperatures, points A, B 
and C (Fig. 6) will move to the right and up. They will reach the p/po = 0.5 line. These 
define the end-points T, p in Fig. 7. We do not continue to higher densities as the simple 
approximation of excluded volume as a means of incorporating interactions between clusters 
becomes progressively worse. If we accept the validity of the simple multifragmentation 
model in the region p/po < 0.5 we will have to conclude that the critical point does not 
exist in the region p/po < 0.5. The same conclusion can be guessed from other published 



work. Multifragmentation with one kind of particles was also studied by Bugaev et al 
This is the same physics problem as considered in [3] but treated in a different mathematical 
framework and these authors considered all densities, not just p/po < 0.5. They found that 
one can identify a critical point at T = Tc = 18.0 MeV, p/po = 1 and pc = oo. At very high 
pressure the model must break down but this is an additional confirmation that the simple 
multifragmentation model in the domain p/po < 0.5 does not contain the critical point. 



VII. ISOTHERMALS IN A TWO-COMPONENT SYSTEM 

Figure 6 gives the isothermals for A^o/^o=l, 1-4 and 1.8 at 6.5 MeV temperature. Drawing 
isothermals for fixed Nq/Zq is physically relevant. We are assuming that we have a very large 
system with given numbers Nq,Zq whose volume can change depending upon the physical 
conditions it is subjected to. If we want to study a different asymmetry we change Nq/Zq 
accordingly and repeat the calculation. To have a complete knowledge, calculations should 
be done for all relevant Nq/Zq. The most asymmetric system we study is Nq/Zq=1.8. 
Of course, since we have no Coulomb force the system with Nq/Zq = a has the same 
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thermodynamic properties as the system with Zq/Nq — a. 

It is however instructive to consider isothermals of two-component systems in a more 
general fashion. In a one-dimensional system there is only one density and an isotherm is 
a line in the p — p plane. Now we have two densities p„ and Pp and isothermals become 
surfaces. Let Pp be the x-axis, p„ the y-axis and p the z-axis, the equation of state at a given 
temperature is a surface in this space. A projection of this surface in two dimensions can 
be made but for a quantitative study it is more convenient to present contours of constant 
p in Pp, pn plane. Such a plot is shown in Fig.8. We consider pressure contours in the region 
bounded by pp — 1.8pn,Pn = l-8pp and (p„ -|- pp)/ po < 0.5. The reasons for choosing these 
boundaries were explained before. 

Roughly speaking, the contours are cither largely radial or circular. Let us first consider 
an uninteresting gas. We assume it consists of only neutrons and protons and unlike in the 
present problem does not form composites. In such a case constant pressure curves would 
be Pn + Pp=constant and these would be straight lines making angle of 7r/4 with the x and y 
axes. Instead we see at low p„ and pp (when one has a gas phase only) not straight lines but 
more like concentric circles. This is because pressure is directly proportional to multiplicity 
(section III) and multiplicity is a function of asymmetry. In our case, composites arc present 
in the gas phase and the number of composites depend upon the asymmetry of the system. 
This causes constant pressure contours in the gas phase to bend from straight lines. We skip 
the details why the lines become like circles. We now try to explain other pressure contours 
which are largely radial. For this, refer back to Fig.6. We have mentioned before that in 
the limit Amax oo the p — p curves would have zero slopes to the right of points A, B and 
C on the isothermals. In such a case constant pressure contour would move exactly radially 
inwards from the boundary p/po — 0.5, would later leave the radial pattern, bend and finish 
at the boundary p„ = l-8pp or pp — 1.8p„ whichever is appropriate. Similar behaviour is 
seen in Fig.8. Thus radial pressure contours reflect regions of co-existence. 

As another example we show in Fig. 9 the pressure contours at T—7.5 MeV. Except near 
the edges of the boundaries, pure gas phase is seen. 
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VIII. ISOSPIN FRACTIONATION 



We illustrate isospin fractionation in multifragmentation model through an example. 
Consider multifragmentation of a neutron rich system: Nq/ Zq=1A and temperature T=6.5 
MeV. At low density the system is in pure gas phase. Following section IV, the gas phase 
consists of light particles with A < 70 and the liquid phase consists of particles with A 
between A^ax-^^ and A^ax- At higher density, both gas and liquid phases are seen (Figs. 
2 and 3). In the present example with Nq/ Zq=1A, we expect that during co-existence the 
neutron to proton ratio in the gas phase will rise above 1.4 and the neutron to proton ratio 
in the liquid phase will fall below 1.4. The reason for this is the symmetry energy which 
preferentially favours formation of larger clusters closest to maximum stability (i.e.,A^ = Z). 
This rise of neutron to proton ratio in the gas phase is illustrated in Fig. 10. Co-existence 
sets in a little beyond p/po=0.1. Till that point is reached the neutron to proton ratio in 
the gas phase is at 1.4, the ratio of the parent system. Then as the density increases the 
ratio increases. 

Fig. 10 also shows that even at very low density the ratio of unbound single neutrons to 
unbound single protons rises very rapidly. But this has got nothing to do with what is called 
isospin fractionation. In fact nothing special happens to this ratio when co-existence sets 
in. It is only if the gas phase is considered to be not just single nucleons but includes light 
particles as well that isospin fractionation becomes an order parameter ii Nq/ Zq ^ 1. 

In the present example, at p/po = 0.35 the neutron to proton ratio in the gas phase is 
1.485. In the liquid phase it is 1.375. 

Isospin fractionation in mean-field theories is treated in 0, [sl. Calculations in the lattice 
gas model can be found in |l^ . 




IX. SUMMARY 

The multifragmentation model, so useful for fitting experimental data in intermediate 
energy collisions, leads naturally to a model of phase transition for nuclear matter. In a 
range of temperature and density first order phase transition occurs. The gas phase and 
the liquid phase can be clearly identified. This is really remarkable. The model of nuclear 
multifragmentation may be unique in this respect. The gas phase consists of light nuclei 
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with A up to about 70. Besides these gas particles, there are large blobs of matter (hquid) 
with mass numbers close to A^ax with A^ax oo. The model is appropriate at subnormal 
nuclear density. Modifications of the simple model are needed to extend the model to higher 
density but this may not be easy. 

Actual nuclear systems as created in heavy ion collisions are finite and in addition have 
Coulomb forces. This makes identification of signals which are finger prints of phase transi- 
tion difficult. This continues to be the subject of intense study and there is large volume of 
literature but this is outside the scope of the present article. 
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FIG. 1: Pressure-density curves for Nq/Zq = 1.4 and T = 6.5 MeV, where the values of Amax used 
are 200, 400, 600, 800 and 1000. Note that in the region of fast rise of pressure with density results 
are insensitive to the value A^ax- In the high density side pressure appears to approach a constant 
value as a function of density as the the value of A^ax is increased. 
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FIG. 2: Plot of A < UA > /Aq as a function of the mass number A for Nq/Zq = 1.0 and T=6.5 
MeV. The sohd hne gives the distribution of composites at p/po = 0.1. There are practically no 
heavy particles, none above ^=70. This is pure gas phase. The dotted line is at p/po = 0.3. Now 
there are both light and heavy {A > 500) particles. This is co-existence. Here and in the rest of 
the figures we used Ajnax=^00. 
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FIG. 3: Same as Fig.2 except that it is for Nq/Zq = 1.8 and the dotted hne is for p/po = 0.38. 
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FIG. 4: Plot of chemical potential as function of density for No/Zq = 1.4 and T = 6.5 MeV. The 
dotted line is the neutron chemical potential /Xn, the dashed line is the proton chemical potential 
Hp and the solid line is = ^/x^ + ^/^p- 
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FIG. 5: Pressure-density isotherms at r=6, 6.5 ,7.0, 7.25 and 7.5 MeV for iVo/Zo = 1.4 and 
A-max = 600. Note that the point of the beginning of co-existence moves up and to the right as the 
temperature increases. 
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FIG. 6: Pressure-density curves at T = 6.5 MeV for three systems with {Nq/Zq) values equal to 1, 
1.4 and 1.8. The points marked A, B and C on the isotherms will give the values of pressure when 
co-exsistence sets in at T = 6.5 MeV for these Nq/Zq values. 
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T (MeV) 

FIG. 7: Phase-coexistence lines in the p — T plane for different values of Nq/Zq. As in Fig. 6 the 
points marked A, B and C gives the value of pressure where co-existence sets in at r=6.5 MeV. 
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p = 0.01 




FIG. 8: Contours of constant pressure p in pn^Pp plane at T = 6.5 MeV for Amax = 600. The 
values of the pressure (in MeV/fm^) are marked against the contours and some are given in the box 
in the upper right corner. The region is bounded by p„ = 1.8pp, pp = 1.8p„ and {pn + Pp)/po < 0.5. 
The line pn = Pp is shown in the middle. 
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FIG. 10: This figure is for r=6.5 MeV and Nq/Zq = 1.4. The left panel shows the rise of the ratio 
of the number of free neutrons to the number of free protons as a function of density. While the rise 
is fast, nothing particularly new happens at the onset of co-cxistcncc. If, however, the gas phase 
is defined to be all particles with ^4 < 70 (this would be consistent with Figs 2 and 3), the ratio of 
neutrons to protons bound in the gas phase remains that of the parent system till co-existence sets 
in (right panel) and then begins to rise. It behaves like an order parameter if the parent system is 
asymmetric. 
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